!     -------------
!     PROGRAM PROG3
!     -------------
!
!     PLANETARY TABLES FROM 1900 TO 2100.
!
!     Ref : JC-GF February 2000.
!
! --- Object -----------------------------------------------------------
!
!     Example for using the binary direct access file 'platab.dat' :
!     - Read the planetary tables.
!     - Computation of equatorial rectangular coordinates of the bodies
!       (geocentic, heliocentric, barycentric frame J2000).
!
!     First Julian date: JD2415020.5.
!     Step length: 29000.5 days.
!     Number of dates: 3.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      character*10 body(11)
      character*12 fich1,fich2
      character*18 frame(3)
!
      dimension r(6)
!
      data body/
     .'Mercury','Venus','E-M Baryc','Mars','Jupiter','Saturn',
     .'Uranus','Neptune','Pluto','Earth','Sun'/
!
      data frame/
     .'Geocentric frame','Heliocentric frame','Barycentric frame'/
!
! --- Parameters -------------------------------------------------------
!
      t0=2415020.5d0
      Step=29000.5d0
      ndat=3
!
! --- Binary direct access file (unit 10) ------------------------------
!
      lu=10
      fich1='platab.dat'
      open (lu,file=fich1,access='direct',recl=800,status='old',
     .      iostat=nerr)
      if (nerr.ne.0) stop ' Error Direct Access File '
!
! --- Results file (unit 20) -------------------------------------------
!
      fich2='prog3.txt'
      open (20,file=fich2)
      write (20,1001)
      write (*,1001)
!
! --- Compute the coordinates ------------------------------------------
!
      do ib=1,11
         t=t0
         ic=mod(ib+2,3)
         if (ic.eq.0) ic=3
         do n=1,ndat
            call RECTCD (t,ib,ic,1,lu,r,ierr)
            if (ierr.ne.0) then
               write (*,1002) ierr
               stop
            endif
            write (20,1003) body(ib),t,frame(ic),r
            write (*,1003)  body(ib),t,frame(ic),r
            t=t+step
         enddo
      enddo
      write (20,1004) fich2
      write (*,1004)  fich2
      stop
!
! --- Formats ----------------------------------------------------------
!
1001  format (/2x,'PLANETARY TABLES : PROG3'/
     .        /2x,'Rectangular coordinates (equator and equinox J2000)')
1002  format (/2x,'ERREUR : ',i2/)
1003  format (/2x,a/2x,'JD:',f16.5,2x,a/
     .       2x,'x :',f16.12,2x,'y :',f16.12,2x,'z :',f16.12,' au'/
     .       2x,'x'':',f16.12,2x,'y'':',f16.12,2x,'z'':',f16.12,' au/d')
1004  format (/2x,'DONE - Results File : ',a/)
!
      end
!
!
!
      subroutine RECTCD (tjd,ibody,icent,ivel,lu,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : JC-GF February 2000.
!
! --- Object -----------------------------------------------------------
!
!     Rectangular coordinates of planets, earth and sun in geocentric,
!     heliocentric or barycentric frame (equator and equinox J2000).
!
! --- Iutput -----------------------------------------------------------
!
!     tjd         Julian date TDB (double real).
!
!     ibody       Body index (integer).
!                 01  : Mercury.
!                 02  : Venus.
!                 03  : Earth-Moon Barycenter.
!                 04  : Mars.
!                 05  : Jupiter.
!                 06  : Saturn.
!                 07  : Uranus.
!                 08  : Neptune.
!                 09  : Pluto.
!                 10  : Earth.
!                 11  : Sun.
!
!     icent       Center index (integer).
!                 1 : Earth        Geocentric coordinates.
!                 2 : Sun          Heliocentric coordinates.
!                 3 : Barycenter   Barycentric coordinates.
!
!     ivel        Velocity index (integer).
!                 0 : without velocities.
!                 1 : with velocities.
!
!     lu          Logical unit of direct access file (integer).
!
! --- Ouput ------------------------------------------------------------
!
!     r(6)        Table of rectangular coordinates (double real)
!                 r(1) : X  position (au).
!                 r(2) : Y  position (au).
!                 r(3) : Z  position (au).
!                 r(4) : X' velocity (au/day).
!                 r(5) : Y' velovity (au/day).
!                 r(6) : Z' velocity (au/day).
!
!     ierr        Error index (integer).
!                 0 :  no error.
!                 10 : file error.
!                 11 : date error.
!                 12 : body index error.
!                 14 : velocity index error.
!
! --- Subroutines ------------------------------------------------------
!
!     PLANET : heliocentric planetary coordinates.
!     EARTH  : geocentric Earth-Moon barycenter coordinates.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension r(6)
      dimension v1(6),v2(6),v3(6),w(6,9),gm(9)
!
! --- Initialization ---------------------------------------------------
!
      data gms/0.2959122082855911D-03/
      data gm/
     .     0.4912547451450812D-10,
     .     0.7243452486162703D-09,
     .     0.8997011374291877D-09,
     .     0.9549535105779258D-10,
     .     0.2825345909524226D-06,
     .     0.8459715185680659D-07,
     .     0.1292024916781969D-07,
     .     0.1524358900784276D-07,
     .     0.2188699765425970D-11/
!
      data tdeb/2415020.0d0/
      data tfin/2487981.0d0/
!
! --- Clear the results tables -----------------------------------------
!
      max=3
      if (ivel.eq.1) max=6
      do i=1,max
         r(i)=0.d0
      enddo
!
! --- Check up the parameters ------------------------------------------
!
      if (ibody.eq.10.and.icent.eq.1.or.
     .    ibody.eq.11.and.icent.eq.2) return
!
      ierr=11
      if (tjd.lt.tdeb.or.tjd.gt.tfin) return
      ierr=12
      if (ibody.lt.1.or.ibody.gt.11) return
      ierr=13
      if (icent.lt.1.or.icent.gt.3) return
      ierr=14
      if (ivel.lt.0.or.ivel.gt.1) return
      ierr=0
!
! --- Select the frame -------------------------------------------------
!
      select case (icent)
!
! --- Geocentric coordinates -------------------------------------------
!
      case (1)
         if (ibody.eq.11) then
            call PLANET (tjd,3,ivel,lu,v1,ierr)
            if (ierr.ne.0) return
            call EARTH (tjd,ivel,v2,ierr)
            if (ierr.ne.0) return
            do i=1,max
               v1(i)=v2(i)-v1(i)
            enddo
         else
            if (ibody.eq.3) then
               call EARTH (tjd,ivel,v1,ierr)
               if (ierr.ne.0) return
            else
               call PLANET (tjd,ibody,ivel,lu,v1,ierr)
               if (ierr.ne.0) return
               call PLANET (tjd,3,ivel,lu,v2,ierr)
               if (ierr.ne.0) return
               call EARTH (tjd,ivel,v3,ierr)
               if (ierr.ne.0) return
               do i=1,max
                  v1(i)=v1(i)-v2(i)+v3(i)
               enddo
            endif
         endif
!
! --- Heliocentric coordinates -----------------------------------------
!
      case (2)
         if (ibody.eq.11) return
         if (ibody.eq.10) then
            call PLANET (tjd,3,ivel,lu,v1,ierr)
            if (ierr.ne.0) return
            call EARTH (tjd,ivel,v2,ierr)
            if (ierr.ne.0) return
            do i=1,max
               v1(i)=v1(i)-v2(i)
            enddo
         else
            call PLANET (tjd,ibody,ivel,lu,v1,ierr)
            if (ierr.ne.0) return
         endif
!
! --- Barycentric coordinates ------------------------------------------
!
      case (3)
         dmas=gms
         do i=1,max
            v1(i)=0.d0
         enddo
         do k=1,9
            call PLANET (tjd,k,ivel,lu,w(1,k),ierr)
            if (ierr.ne.0) return
            dmas=dmas+gm(k)
            do i=1,max
               v1(i)=v1(i)-gm(k)*w(i,k)
            enddo
         enddo
         do i=1,max
            v1(i)=v1(i)/dmas
         enddo
         if (ibody.ne.11) then
            if (ibody.ne.10) then
               do i=1,max
                  v1(i)=v1(i)+w(i,ibody)
               enddo
            else
               call EARTH (tjd,ivel,v2,ierr)
               if (ierr.ne.0) return
               do i=1,max
                  v1(i)=v1(i)+w(i,3)-v2(i)
               enddo
            endif
         endif
!
      endselect
!
! --- Results ----------------------------------------------------------
!
      do i=1,max
         r(i)=v1(i)
      enddo
!
      return
      end
!
!
!
      subroutine PLANET (tjd,ipla,ivel,lu,r,ierr)
!-----------------------------------------------------------------------
!
!     PLANETARY TABLES FROM 1900 TO 2100.
!
!     Ref : JC-GF February 2000.
!
! --- Object -----------------------------------------------------------
!
!     Equatorial rectangular coordinates J2000 of the planets computed
!     with the planetary tables stored in a binary direct access file.
!
! --- Input ------------------------------------------------------------
!
!     tjd         Julian date TDB  (double real).
!                 Time span of validity : JD2415020.5 - JD2487981.5.
!                 tdeb = 2415020.0    31/12/1899 12h.
!                 tfin = 2487981.0    04/10/2099 12h.
!
!     ipla        Planet index (integer).
!                 1 : Mercury.
!                 2 : Venus.
!                 3 : Earth-Moon Barycenter.
!                 4 : Mars.
!                 5 : Jupiter.
!                 6 : Saturn.
!                 7 : Uranus.
!                 8 : Neptune.
!                 9 : Pluto.
!
!     ivel        Velocity index (integer).
!                 0 : without velocities.
!                 1 : with velocities.
!
!     lu          Logical unit of direct access file (integer).
!
! --- Output -----------------------------------------------------------
!
!     r(6)        Table of rectangular coordinates (double real).
!                 r(1) : X  position (au).
!                 r(2) : Y  position (au).
!                 r(3) : Z  position (au).
!                 r(4) : X' velocity (au/day).
!                 r(5) : Y' velovity (au/day).
!                 r(6) : Z' velocity (au/day).
!
!     ierr        Error index (integer).
!                 0 :  no error.
!                 10 : file error.
!                 11 : date error.
!                 12 : planet index error.
!                 14 : velocity index error.
!
! --- Binary direct access file ----------------------------------------
!
!     The planetary tables are stored in a binary direct access file.
!     This file has to be opened in the monitor (1 record = 800 bytes).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension r(6)
      dimension v(6,9)
      dimension fq(1500),cf(9108)
      dimension tzero(9),kblk(9),dt(9),mx(9),nf(9,0:2)
      dimension lfq(9),lcf(9),nr0(9),nrf(9),nrt(9)
!
! --- Initialization ---------------------------------------------------
!
      data tdeb/2415020.0d0/
      data tfin/2487981.0d0/
      data tzero/7*2415020.5d0,2396758.5d0,2341972.5d0/
      data kblk/20,8,5,5,1,1,1,1,1/
      data dt/3650.d0,9120.d0,2*14630.d0,3*73080.d0,91400.d0,146120.d0/
      data mx/2,2,2,2,2,2,2,2,2/
      data nf/224,178,208,183,134,117,087,080,079,
     .        013,021,025,024,024,022,018,018,022,
     .        002,001,002,001,004,002,001,004,006/
      data lfq/1,240,440,675,883,1045,1186,1292,1394/
      data lcf/1,1447,2659,4081,5341,6325,7183,7831,8455/
      data nr0/1,304,410,488,556,568,579,588,597/
      data nrf/3,2,3,3,2,2,2,2,2/
      data nrt/15,13,15,13,10,9,7,7,7/
!
      max=3
      if (ivel.eq.1) max=6
      do i=1,max
         r(i)=0.d0
      enddo
!
! --- Check up the parameters ------------------------------------------
!
      ierr=11
      if (tjd.lt.tdeb.or.tjd.gt.tfin) return
      ierr=12
      if (ipla.lt.1.or.ipla.gt.9) return
      ierr=14
      if (ivel.lt.0.or.ivel.gt.1) return
      ierr=0
!
! --- Read the series --------------------------------------------------
!
      nblk=(tjd-tzero(ipla))/dt(ipla)+1
      if (nblk.lt.1) nblk=1
      if (nblk.gt.kblk(ipla)) nblk=nblk-1
!
      mf=nf(ipla,0)+nf(ipla,1)+nf(ipla,2)
      mt=(mf*2+4)*3
      lf=lfq(ipla)-1
      lc=lcf(ipla)-1
!
      nrec=nr0(ipla)-1
      do n=1,nrf(ipla)
         nrec=nrec+1
         jmin=lf+(n-1)*100+1
         jmax=min0(lf+mf,jmin+99)
         read (lu,rec=nrec,err=100) (fq(j),j=jmin,jmax)
      enddo
!
      nrec=nrec+(nblk-1)*nrt(ipla)
      do n=1,nrt(ipla)
         nrec=nrec+1
         jmin=lc+(n-1)*100+1
         jmax=min0(lc+mt,jmin+99)
         read (lu,rec=nrec,err=100) (cf(j),j=jmin,jmax)
      enddo
!
! --- Compute the coordinates ------------------------------------------
!
      t0=tzero(ipla)+(nblk-1)*dt(ipla)
      x=2.d0*(tjd-t0)/dt(ipla)-1.d0
      fx=dt(ipla)*x/2.d0
      kc=lcf(ipla)-1
!
! --- Compute the positions --------------------------------------------
!
      do iv=1,3
!
         v(iv,ipla)=0.d0
         wx=1.d0
!
         do i=0,3
            kc=kc+1
            v(iv,ipla)=v(iv,ipla)+cf(kc)*wx
            wx=wx*x
         enddo
!
         wx=1.d0
         kf=lfq(ipla)-1
!
         do m=0,mx(ipla)
            ws=0.d0
            do i=1,nf(ipla,m)
               kf=kf+1
               f=fq(kf)*fx
               c=cos(f)
               s=sin(f)
               ws=ws+cf(kc+1)*c+cf(kc+2)*s
               kc=kc+2
            enddo
            v(iv,ipla)=v(iv,ipla)+ws*wx
            wx=wx*x
         enddo
!
      enddo
!
      do iv=1,3
         v(iv,ipla)=v(iv,ipla)/1.d10
         r(iv)=v(iv,ipla)
      enddo
      if (ivel.eq.0) return
!
! --- Compute the velocities -------------------------------------------
!
      wt=2.d0/dt(ipla)
      kc=lcf(ipla)-1
!
      do iv=4,6
!
         v(iv,ipla)=0.d0
         wx=1.d0
         kc=kc+1
!
         do i=1,3
            kc=kc+1
            v(iv,ipla)=v(iv,ipla)+i*cf(kc)*wx
            wx=wx*x
         enddo
         v(iv,ipla)=wt*v(iv,ipla)
!
         wx=1.d0
         kf=lfq(ipla)-1
!
         do m=0,mx(ipla)
            do i=1,nf(ipla,m)
               kf=kf+1
               fw=fq(kf)
               f=fw*fx
               c=cos(f)
               s=sin(f)
               ctw=cf(kc+1)
               stw=cf(kc+2)
               v(iv,ipla)=v(iv,ipla)+fw*(stw*c-ctw*s)*wx
               if (m.gt.0)
     .         v(iv,ipla)=v(iv,ipla)+m*wt*(ctw*c+stw*s)*wy
               kc=kc+2
            enddo
            wy=wx
            wx=wx*x
         enddo
!
      enddo
!
      do iv=4,6
         v(iv,ipla)=v(iv,ipla)/1.d10
         r(iv)=v(iv,ipla)
      enddo
!
      return
!
! --- File error -------------------------------------------------------
!
100   continue
!
      ierr=10
      do iv=1,max
         r(iv)=0.d0
      enddo
!
      return
      end
!
!
!
      subroutine EARTH (tjd,ivel,r,ierr)
! ----------------------------------------------------------------------
!
!     Ref : JC-GF February 2000.
!
! --- Object -----------------------------------------------------------
!
!     Geocentric rectangular coordinates of Earth-Moon barycenter.
!     Frame: mean equator and equinox J2000.
!     Source: ELP2000-82 (truncated series).
!
! --- Input ------------------------------------------------------------
!
!     tjd :       Julian date TDB (double real).
!
!     ivel        Velocity index (integer).
!                 0 : without velocities.
!                 1 : with velocities.
!
! --- Output -----------------------------------------------------------
!
!     r(6) :      Table of rectangular coordinates (double real).
!                 r(1) : x position (ua).
!                 r(2) : y position (ua).
!                 r(3) : z position (ua).
!                 r(4) : x vitesse  (ua/j).
!                 r(5) : y vitesse  (ua/j).
!                 r(6) : z vitesse  (ua/j).
!
!     ierr        Error index (integer).
!                 0 :  no error.
!                 14 : velocity index error.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension r(6)
      dimension v(6),n1(3),n2(3)
      dimension c(138),s(138),f(138)
!
! --- Initialization ---------------------------------------------------
!
      data n1/01,44,93/
      data n2/43,92,138/
!
! --- Coefficients of cosines (10*10 au) -------------------------------
!
      data c/
     .-244075.,  -2965.,   8528.,   2345.,  -2486.,   1426.,    527.,
     .    -43.,   -393.,    394.,   -218.,     73.,     91.,   -173.,
     .     25.,    -20.,     75.,     72.,      6.,     72.,    -40.,
     .    -58.,     56.,    -53.,     46.,    -44.,     -5.,      0.,
     .     -1.,    -12.,     21.,     -4.,     -4.,     -2.,      9.,
     .      8.,     10.,      2.,    -10.,    -12.,    -11.,     10.,
     .    -10.,-176962., -23344., -11109.,   -922.,  -4118.,    714.,
     .  -1135.,   -601.,    299.,    564.,   -311.,    261.,   -251.,
     .    254.,    229.,    213.,   -179.,     57.,    -19.,   -125.,
     .   -113.,    -87.,    -52.,     75.,     16.,     -5.,    -42.,
     .     -4.,    -10.,      8.,    -19.,      9.,     40.,     29.,
     .    -25.,     11.,     19.,    -12.,     -5.,     18.,    -15.,
     .    -16.,     13.,     13.,      9.,     -1.,     -3.,     -5.,
     .     -1., -76714.,  25611., -10120.,   -400.,   1387.,  -1785.,
     .    310.,    580.,   -492.,   -527.,     44.,    130.,    244.,
     .   -135.,    113.,    -38.,    110.,     92.,    -78.,     25.,
     .    -54.,    -26.,    -49.,    -38.,     27.,    -23.,     32.,
     .      2.,     -2.,      1.,    -18.,     -2.,    -23.,     -4.,
     .      3.,     -8.,      4.,    -11.,     17.,      3.,    -11.,
     .     13.,      8.,    -11.,     -7.,      4./
!
! --- Coefficients of sines (10*10 au) ---------------------------------
!
      data s/
     . 192874.,  25444.,   1005.,   4489.,   -778.,   1238.,   -326.,
     .   -614.,    339.,   -285.,   -276.,   -232.,    195.,    -63.,
     .    136.,    124.,     95.,     57.,    -82.,      5.,     45.,
     .      4.,     11.,     -9.,     20.,    -10.,    -44.,    -32.,
     .     27.,    -20.,      6.,    -20.,     17.,     17.,    -14.,
     .    -14.,      4.,    -11.,    -10.,      3.,      6.,     -7.,
     .      4.,-223938.,  -2720.,    635.,   7824.,   2151.,  -2281.,
     .   1309.,   -675.,    483.,    -40.,   -360.,    362.,    327.,
     .   -200.,    204.,     67.,     84.,   -159.,   -145.,     23.,
     .    -18.,     69.,     66.,      5.,     65.,     66.,    -36.,
     .    -53.,     51.,    -48.,     42.,    -40.,     -5.,      0.,
     .     -1.,    -19.,    -11.,     17.,     20.,     -4.,     -4.,
     .     -2.,      9.,      7.,     -9.,    -13.,    -11.,    -10.,
     .     11., -97079.,  -1464.,  -1179.,   3392.,   1557.,    933.,
     .   -989.,   -754.,    567.,   -470.,    334.,    210.,    -17.,
     .   -156.,    157.,   -151.,    -87.,     29.,     36.,    -69.,
     .     10.,     43.,     -8.,     30.,    -39.,     29.,      2.,
     .     29.,     29.,    -25.,    -16.,    -23.,      3.,     22.,
     .    -21.,     18.,    -17.,     13.,     -2.,     14.,     -8.,
     .      0.,     -9.,      0.,     -8.,      9./
!
! --- Frequencies (rad/d) ----------------------------------------------
!
      data f/
     . 0.2299708345453799, 0.0019436907548255, 0.4579979783362081,
     . 0.0324605575663244,-0.1955665862245038, 0.4274811115247091,
     .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655,
     . 0.2471728045705681, 0.6860251221267625, 0.2604877013568789,
     . 0.0496625275912389,-0.1783646161993155, 0.0172021241604381,
     . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209,
     .-0.0152582792703628, 0.4407960083110198, 0.8835353991060917,
     . 0.1994539677341547, 0.2662248529612594, 0.4751999483613963,
     . 0.4427395449305955,-0.0037934608495551, 0.6383062852903491,
     . 0.4637351299405887, 0.1937168161297741, 0.0152585875411362,
     .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377,
     . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383,
     . 0.2300338378809035,-0.2299078312098563, 0.4446830815498973,
     . 0.8530185322948665, 0.4885148451477070, 0.0381977091707050,
     . 0.2147124011397673, 0.2299708345453799, 0.0019436907548255,
     . 0.2308957195928816, 0.4579979783362081, 0.0324605575663244,
     .-0.1955665862245038, 0.4274811115247091,-0.0028685758023272,
     .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655,
     . 0.2471728045705681, 0.1946417011770021, 0.6860251221267625,
     . 0.4589228633837098, 0.2604877013568789, 0.0496625275912389,
     .-0.1783646161993155,-0.0333854426135524, 0.0172021241604381,
     . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209,
     .-0.0152582792703628, 0.4284059965722108, 0.4407960083110198,
     . 0.8835353991060917, 0.1994539677341547, 0.2662248529612594,
     . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551,
     . 0.6383062852903491, 0.4637351299405887, 0.1937168161297741,
     . 0.6564331403627652, 0.0152585875411362, 0.1774397311520876,
     .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377,
     . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383,
     . 0.4446830815498973, 0.6869500071742641, 0.8530185322948665,
     . 0.4885148451477070, 0.2251585679885010, 0.2299708345453799,
     . 0.2308957195928816, 0.0019436907548255, 0.4579979783362081,
     .-0.0028685758023272, 0.0324605575663244,-0.1955665862245038,
     . 0.1946417011770021, 0.4274811115247091, 0.4589228633837098,
     .-0.0333854426135524,-0.2318206046403833, 0.6555082553155372,
     . 0.2127688645204655, 0.2471728045705681, 0.4284059965722108,
     . 0.6860251221267625, 0.2604877013568789, 0.0496625275912389,
     .-0.1783646161993155, 0.0172021241604381, 0.6564331403627652,
     . 0.0191456607800137,-0.2260834530360027, 0.1774397311520876,
     . 0.4102791414995209,-0.0152582792703628, 0.6869500071742641,
     . 0.4407960083110198, 0.2251585679885010, 0.8835353991060917,
     . 0.1994539677341547, 0.4226688449678302, 0.2662248529612594,
     . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551,
     . 0.2118436712021903, 0.6383062852903491,-0.0505874126387406,
     .-0.2614125864043805, 0.4637351299405887, 0.0200705458272416,
     . 0.1937168161297741, 0.0143333942228611,-0.0181270092079398/
!
! --- Clear the results tables -----------------------------------------
!
      max=3
      if (ivel.eq.1) max=6
      do i=1,max
         r(i)=0.d0
      enddo
!
! --- Check up the velocity index --------------------------------------
!
      ierr=14
      if (ivel.lt.0.or.ivel.gt.1) return
!
! --- Compute the coordinates ------------------------------------------
!
      t=tjd-2451545.0d0
!
      do iv=1,3
         v(iv)=0.d0
         if (ivel.eq.1) v(iv+3)=0.d0
         do n=n1(iv),n2(iv)
            x=f(n)*t
            cx=cos(x)
            sx=sin(x)
            v(iv)=v(iv)+c(n)*cx+s(n)*sx
            if (ivel.eq.1) v(iv+3)=v(iv+3)-c(n)*f(n)*sx+s(n)*f(n)*cx
         enddo
         v(iv)=v(iv)/1.d10
         r(iv)=v(iv)
         if (ivel.eq.1) then
            v(iv+3)=v(iv+3)/1.d10
            r(iv+3)=v(iv+3)
         endif
      enddo
!
      ierr=0
      return
      end
